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DEBLURRING  GAUSSIAN  BLUR 


Robert  A.  Hummel 
B.  Kimia 
Steven  W.  Zucker 

Abstract 

Gaussian  blur,  or  convolution  against  a  Gaussian  kernel,  is  a  common  model 
for  image  and  signal  degradation.  In  general,  the  process  of  reversing 
Gaussian  blur  is  unstable,  and  cannot  be  represented  as  a  convolution  filter 
in  the  spatial  domain.  If  we  restrict  the  space  of  allowable  functions  to  poly¬ 
nomials  of  fixed  finite  degree,  then  a  convolution  inverse  does  exist.  We 
give  constructive  formulas  for  the  deblurring  kernels  in  terms  of  Hermite 
polynomials,  and  observe  that  their  use  yields  optimal  approximate  deblur¬ 
ring  solutions  among  the  space  of  bounded  degree  polynomials.  The  more 
common  methods  of  achieving  stable  approximate  deblurring  include  restric¬ 
tions  to  band-limited  functions  or  functions  of  bounded  norm. 


1.  Introduction 

Given  an  image  or  signal,  the  realization  of  any  system  for  processing  it  must 
introduce  some  amount  of  degradation.  Since  there  may  be  several  stages  each  con¬ 
tributing  to  the  degradation,  the  composition  is  often  modeled  as  a  Gaussian  blur¬ 
ring  operation.  We  consider  spatially  invariant  Gaussian  convolution  defined  as  fol¬ 
lows.  For  a  bounded  measurable  input  function  / (x)  defined  for  xeIR",  then  the 
observed  blurred  output  is  given  by 

h(x)  =  A'Or- {,»)/( €)</$, 

where 


K(x,t) 


_ 1 _ e 

(4170"  2 


and  f  is  a  fixed  positive  value  parameterizing  the  extent  of  the  blur, 
estimate  /(.v)  when  only  h(x)  and  the  amount  of  blur  t  are  known. 


We  wish  to 


It  would  be  especially  nice  to  formulate  deblurring  as  a  convolution  operation, 
so  that 


/  =  D(-,0  *  h. 

In  general,  a  universal  deblurrinc  kernel  Dt  v,n  does  not  exist  Howexer,  il  sutti- 
cicnt  restrictions  are  placed  on  the  domain  of  permissible  (unctions  /,  then  deblur¬ 
ring  kernels  can  exist 

Our  interest  in  deblurring  is  motivated  by  two  concerns.  First,  deblurring  is  o! 
significant  practical  importance  in  mans  image  processing  systems,  lur  example  it-. 
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computerized  tomography  [1].  There  are  also  applications  in  physiological  optics, 
such  as  the  de-focusing  that  automatically  takes  place  for  objects  outside  the  depth 
of  field  of  an  accommodated  eye. 

A  second  motivation  is  provided  by  a  desire  to  study  stability  of  image 
representations.  The  mathematical  analysis  of  an  image  representation  must  include 
a  study  of  the  continuity  and  stability  of  the  transformation.  Reconstruction  methods 
are  particularly  useful  for  studying  the  stability.  While  blurred  versions  of  an  origi¬ 
nal  signal  form  a  classical  unstable  representation,  many  intermediate-level  transfor¬ 
mations  of  image  data  nonetheless  involve  some  degree  of  filtering  by  blurring.  For 
example,  representations  involving  zero-crossings  of  Gaussian-filtered  Laplacians  of 
images  [2],  as  well  as  many  other  pyramid  schemes  [3, 4, 5, 6],  involve  Gaussian 
blur.  Instabilities  in  the  representation  may  not  be  important  if  approximate  or 
pseudoinverse  reconstruction  methods  (see,  e.g.,  [7])  can  be  found  that  make  expli¬ 
cit  the  assumptions  concerning  the  input  data.  In  this  paper,  we  present  an 
approximate-inverse  method  involving  polynomial  approximations. 

2.  The  Heat  Equation 

2.1.  Diffusion 

There  is  a  fundamental  connection  between  Gaussian  blurring  and  the  heat 
equation.  Consider  a  rod  of  infinite  length  onto  which  an  impulse  of  heat  is  placed. 
As  time  evolves,  the  heat  will  diffuse  and  the  original  impulse  will  spread  out.  By 
elementary  physics,  the  resulting  distribution  will  approximate  a  Gaussian  whose 
width  depends  on  the  elapse  time  (  see,  e.g.,  the  Feynman  lectures,  [8]).  By  super¬ 
position,  the  model  for  the  temperature  distribution  along  the  rod  at  any  given  time 
is  the  initial  temperature  distribution  convolved  with  a  Gaussian.  The  diffusion  pro¬ 
cess  effectively  convolves  the  initial  distribution  by  a  Gaussian  whose  spread 
depends  on  how  much  time  has  evolved.  This  is  the  physically  realized  solution  to 
the  heat  equation,  which  can  be  formulated  as  follows  (  [9]).  Given  fix)  piecewise 
continuous  and  bounded,  find  h(x,t)  bounded  and  C2  for  r>0  satisfying 

-^-(x,f)  =  A/j(x,r),  xcIR",  r > 0; 

dt 

h(x,t)-f(x0 )  as  (.r,0-Uo.O),  xocIR",  t> 0. 

We  denote  the  operator  that  takes  /  to  h(,t)  by  flf,  i.e.,  h(x,t)  =  (fl, fH  v)  The 
solution  is  given  by 

(fl,f)(x)  =  JrRnK(y-x,l)f(y)dy, 

where  K  is  as  defined  before.  When  restricted  to  a  Hilbert  space  such  ^ 
/.2(lRn),  li,  becomes  a  symmetric  bounded  linear  operator.  We  will  generally  inter¬ 
pret  fix)  as  an  image. 

The  space  of  functions  that  can  be  blurred  in  this  way  is  very  large  Indeed, 
the  condition  that  fix)  is  bounded  can  be  weakened.  The  solution  is  still  given  r-v 
convolution  against  the  “source  kerne!'’  K(x,t),  [10],  hi  j )  =  K(  J)*f  The  source 
kernel  is  the  fundamental  solution  to  the  heat  equation  on  the  unbounded  d  ••  . 

IR  n  We  also  note  that  the  blurring  operator  satisfies  a  semigroup  property  r:  .e 
ll  a  (unction  fix)  has  been  blurred  for  some  time  /•  by  II.  .  and  the  rev.i!:. 
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function  is  blurred  by  ft,2,  the  end  result  is  the  same  as  blurring  f(x)  for  a  time 
r i  +  •  That  is,  ft, oft,  =  ft,+/.  The  two  Gaussian  blurring  operators,  each  of  which 

may  have  its  own  physical  justification,  results  in  one  composite  Gaussian  operator. 
Indeed,  by  the  central  limit  theorem,  other  blurring  operators  will  also  compose  into 
approximate  Gaussians  when  iterated. 

2.2.  Deblurring 

Deblurring  is  the  inverse  problem  to  blurring,  and  can  be  modeled  as  a  diffu¬ 
sion  process  running  backwards  in  time.  Formally,  the  problem  of  reconstructing 
/(*)  given  a  blurred  function  h(x)  and  a  blurring  amount  t> 0  is  the  inverse  heat 
equation  problem,  and  poses  technical  difficulties  not  present  in  the  forward  heat 
equation  problem. 

First,  finding  an  inverse  to  ft,  presupposes  that  ft,  is  one-to-one.  In  fact,  the 
blurring  operator  is  one-to-one  providing  minor  restrictions  are  placed  on  the 
domain  of  ft,.  However,  without  certain  growth  restrictions,  it  is  possible  to  find 
distinct  functions /and /satisfying  ft,/ =  ft,/,  (see  [11]).  Second,  a  solution /to  the 
problem  ft,/  =  h,  given  h,  exists  only  if  h  is  sufficiently  smooth.  In  general,  an 
inverse  can  not  be  found,  and  even  if  h  is  sufficiently  smooth,  an  arbitrarily  small 
change  can  destroy  the  smoothness.  John  [12]  discusses  the  technical  conditions 
needed  for  the  existence  of  an  inverse.  Finally,  in  a  general  function  space  the 
deblurring  problem_is  horribly  ill-conditioned.  This  means  that  there  can  exist  pairs 
of  functions  /  and  /  that  are  arbitrarily  far  apart  whose  image_s  under  ft,  are  arbi¬ 
trarily  close.  The  prototvpic  example  is  f(x)  =  /4sin(u)x)  and  /(*)=  0  in  one  space 
dimension.  Then  (ft,/)(.v)  =  Ae  "“‘'sin(  wr) ,  which  for  u>  large  can  be  very  close  to 

ft,/  =  0. 

Deblurring  can  be  understood  somewhat  better  in  terms  of  the  Fourier 
transform.  If  we  denote  the  Fourier  transform  of  a  function  g(x)  by  g(i o),  then  the 
blurring  operator  ft,  is  a  multiplier  operator  given  by 

(ft, /)*(<*>)  =  <?  -ulf/(<*>) . 

By  means  of  this  formula,  ft,  can  be  extended  to  operate  on  the  class  of  temperate 
distributions  S'  of  Fourier  transformable  distributions  [13].  In  particular,  ft,/  is 
defined  for  any  polynomial  /.  Further,  the  formula  shows  that  ft,  is  one-to-one  on 
any  class  of  Fourier  transformable  functions.  Moreover,  our  earlier  observation 
that  deblurring  can  not  generally  be  represented  by  a  convolution  kernel  can  be 
observed  from  the  formula,  since  although 

/(w)  = 

a  general  convolution  formula  is  not  possible  since  eu‘r  is  not  the  Fourier  transform 
of  any  tempered  distribution 

These  difficulties  would  tend  to  make  one  pessimistic  about  accomplishing 
image  deblurring.  and  in  particular  about  discovering  deblurring  kernels  However, 
deblurring  is  a  common  operation,  and  is  typically  accomplished  b\  giving  the  prob¬ 
lem  a  variational  formulation,  which  can  lead  to  a  well  conditioned  problem  We 
describe  several  variational  formulations  in  the  next  section,  and  present  deblurring 
kernels  for  polynomial  domains  in  the  subsequent  section 
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3.  Variational  Formulation 

We  choose  a  normed  space  M  and  a  closed  convex  subset  X  £  X  so  that  ft, 
may  be  regarded  as  an  operator  Cl,  :  X  -  X,  for  all  0;  (note  that  Clo  is  simply  the 
identity  operator.)  We  may  then  pose  the  deblurring  problem  in  the  following 
form: 

Given  A€X,  t> 0,  find  /£X  minimizing  \\Cl,f  —  h\\. 

If  ft,  is  one-to-one  and  onto  over  X ,  then  the  solution  /  given  h  is  precisely  the 
inverse  image  of  h  under  ft,,  so  that  the  minimization  gives  a  zero  norm.  In  gen¬ 
eral,  however,  X  is  restricted  in  such  a  way  that  ft,  maps  into  X,  and  so  the  solu¬ 
tion  /is  a  pseudoinverse  of  h. 

Difficulties  arise  because  the  operator  ft,  on  the  domain  X  is  in  general  ill- 

conditioned.  There  are  various  approaches  that  one  can  take  to  find  a  good 

deblurred  signal  /  stably  from  a  given  h.  A  standard  approach  is  to  restrict  X  to  a 

sufficiently  small  set.  We  mention  five  possibilities,  and  our  subsequent  treatment 

will  use  one  of  these  approaches  (the  second  one)  as  a  point  of  departure. 

* 

(1)  X  =  L2 ,  and  X  =  {/€  X  | /(to)  =  0  for  |  u>  |  >  p},  for  some  fixed  constant  p. 
Restricting  to  the  space  of  band-limited  signals  (with  a  specified  cut-off)  allows 
stable  deblurring  of  a  blurred  signal  h  by  means  of  the  formula 

a  (e “~'4/i  ( u>)  |u>  |<p 

|u>  |>p’ 

The  function  / can  also  be  writtenas  a  convolution  against  h:  f  -  Afp*/i,  where 
the  Fourier  transform  of  A'p  is  e for  |u>  |<p,  and  zero  elsewhere.  This  is  a 
standard  method  for  deblurring,  although  it  is  well  known  that  A'p  “rings”  over 
a  large  spatial  extent. 

In  a  discrete  setting,  a  similar  (discrete)  convolution  deblurring  kernel  can  be 
formulated,  yielding  appropriate  band-limited  discrete  approximations  to  origi¬ 
nal  signals.  It  turns  out  that  this  method  is  completely  equivalent  to  computing 
a  pseudoinverse  of  the  matrix  representing  the  blurring  operator  by  means  of  a 
singular  value  decomposition. 

(2)  X  =  Ll(e  X'dx),  and  X  =  Polynomials  of  degree  S'  or  less  (fixed  S’),  which 
we  will  denote  by  Pfj-  We  will  see  (in  the  next  section)  that  ft,  is  closed  on  P\ 
Thus  the  pscudoinverse  problem  can  be  solved  as  follows  First,  h  is  projected 
onto  P/v  by  the  linear  orthogonal  projection  operator  of  X  into  Ps  The  result¬ 
ing  polynomial  is  deblurred  by  the  convolution  kernel  D\  to  be  defined  in  Sec¬ 
tion  4  to  yield  the  solution  polynomial  /. 

The  problem  with  this  possibility  is  that  images  arc  typically  far  more  genera! 
than  polynomials  of  degree  S'  Thus  to  get  even  moderate!)  reasonable  deblur¬ 
ring,  S'  has  to  be  very  large,  and  then  applying  D\  becomes  difficult  and 
numerically  unstable.  Although  we  will  make  use  of  the  deblurring  kernels  ()., 
designed  to  deblur  polynomials,  our  implementations  (section  5)  are  more  gen¬ 
eral  than  finding  the  pscudoinverse  / among  the  class  P\ 

(  1 1  A  //’(IHl.andX  =  {/(A'|/?~0}  This  situation,  studied  by  John  1 12]  lead- 


vnmfi 
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"blurry."  Specifically,  suppose  T>t,  and  h  m  Clrg,  some  g iM.  Then  the  prob¬ 
lem  is  to  find  f£M  such  that  ft,/  **  h.  John  studies  bounds  on  the  deblurring 
error,  where  deblurring  is  accomplished  by  exactly  the  same  linear  process,  to 
be  discussed  in  in  the  next  section.  That  is,  he  constructs  an  approximation  / to 
/by  convolving  h  with  a  scaled  version  of  the  kernel  Dy  given  in  Section  4.  The 
result  is  that  the  error  in  reconstructing  /,  \\f  -  f\ |,  can  be  controlled  to  depend 
''/'ntinuously  on  the  error  in  representing  h.  That  is,  small  errors  in  represent¬ 
ing  h  can  lead  to  errors  in  representing  /,  but  the  maximum  size  of  the  errors 
can  be  bounded.  Interestingly,  unlike  customary  notions  of  stability,  the 
dependence  is  not  linear.1 

(4)  J!/  =  L2(1R),  and  W  =  {/|  Os/sA/},  for  some  fixed  M.  With  M  *  ®,  this  is 
nearly  the  same  as  case  (2).  Now,  however,  we  consider  the  possibility  of  non¬ 
linear  deblurring  methods.  Peleg  [14]  has  implemented  a  deblurring  scheme 
base  on  a  conjugate  gradient  iterative  minimization  of  j|ft,/  -  h\\,  constrained 
by  f£M.  The  constraints  are  handled,  in  Peleg’s  case,  by  remapping  the  inter¬ 
val  [0..M]  to  [  —  ao, oo],  and  then  solving  an  unconstrained  minimization  prob¬ 
lem.  By  limiting  the  number  of  iterations,  they  obtain  only  an  approximate 
solution,  although  the  results  look  very  good.  They  don’t  study  the  stability 
question,  but  one  would  expect  the  same  kind  of  nonlinear  stability  for  partial 
deblurring  as  discovered  by  John. 

(5)  M  =  L2(IR),  M  =  {/€P2|  |l/j|5A/},  some  fixed  M.  This  case  is  treated  by 
Carasso  et.  al.  [15].  They  give  a  relatively  simple  nonlinear  deblurring 
method,  making  use  of  Fourier  transforms,  to  solve  the  variational  problem. 
The  method  is  not  iterative.  They  also  study  the  stability,  and  obtain  the  same 
kind  of  stability  estimates  (for  partial  deblurring)  as  John. 

An  alternative  approach  to  obtaining  stable  deblurring  is  to  begin  by  specifying 
the  algorithmic  form  of  the  deblurring  method,  and  to  optimize  with  respect  to  a 
statistical  norm.  For  example,  we  can  insist  on  a  convolution  kernel  for  deblurring, 
and  seek  a  kernel  k  minimizing 

E{\\k*aj  -  /1|}, 

where  £{  }  is  an  expectation  operator  which  presupposes  some  distribution  of 
unblurred  functions  /.  Other  operators,  such  as  worst-case  norm,  are  also  possible. 
Such  methods  are  studied  in  the  province  of  information-based  complexity[16]. 

If  the  distribution  of  /’ s  is  concentrated  on,  or  limited  to,  polynomials  in  Py , 
then  all  the  functions  ft,/  will  also  be  polynomials,  and  we  expect  that  the  optimal 
deblurring  kernel  k  will  be  the  one  that  deblurs  polynomials  in  Py  (namely,  £>y, 
given  in  the  Theorem  of  Section  4).  If  the  distribution  consists  of  functions  that  are 
well-approximated  locally  by  polynomials  in  Py.  the  optimal  kernel  won’t  change 
much. 

In  the  next  section,  we  present  the  kernel  Dy  such  that  Dy*ft,/  =  /  for  all 
/ € Py ,  (with  t  =  1/4).  Since  ft,  is  closed  on  Py,  this  kernel  may  be  used  to  deblur 


1  Suppose  that  (i  is  the  representation  of  h,  and  that  f  it  the  appronmate  reconstruction  of  /  using  h  Usual 
notions  of  stability  would  require  \'J-f  "sj  -Vh-h  ,  using  appropriate  (and  perhaps  different)  norms  The  non¬ 
linear  stability  that  is  used  in  this  case,  however,  asserts  that  | j-f  i  h-h'  •  *,  for  some  integer  Thus 

in  order  to  achieve  an  accuracy  of  «  in  the  reconstruction  of  /.  it  is  necessary  to  represent  h  to  an  accuracy  of 
(1  C  )«*  This  might  be  viewed  as  polynomial  stability 
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all  polynomials  h^Ps-  However,  motivated  by  John’s  studies  (case  (2)  above)  and 
the  comments  on  optimal  kernels,  we  will  apply  Ds  to  functions  h  that  are  not  in 
fact  polynomials. 

The  results  may  be  analyzed  in  terms  of  local  approximations  by  polynomials. 
Suppose  that  the  initial  unblurred  data  /  is  written 

f=P  +  «. 

where  piPs  and  e  is  an  error  term.  Specifically,  p  should  be  the  projection  of  / 
onto  Pn  in  the  space  L2(e~ *  dx),  so  that  p  is  a  good  polynomial  approximation  to  / 
near  the  origin,  but  may  differ  from  / significantly  away  from  the  origin.  Applying 
the  deblurring  kernel  Ds r  to  the  blurred  version  of  /  yields  the  approximate  deblur¬ 
ring 


Thus 


/  =  P  +  DN*£ltf.. 


/-/“«”  DN*D,€, 


whose  norm  (in  L  (e  X‘dx))  can  be  bounded  by  C-||e||.  Since  the  norm  measures 
errors  only  locally,  the  result  is  that  we  have  stable,  accurate  deblurring  for  signals 
that  are  well-approximated  locally  by  polynomial  data. 

4.  Polynomial  Domains 

The  monomials  {1,jc,  •  •  •  ,xN}  form  a  basis  for  PN.  If  this  basis  is  orthonor- 
malized  with  respect  to  the  inner  product 

(J.g)  =  J  f(x)g(x)e~x‘dx, 

then  the  basis  of  Hermite  polynomials  {Ho, Hi,  result.  The  Hermites  can  be 

represented  explicitly: 

ln/2j  (2r)n~2m 

ffn(x)  =  m  2  (~Dm  \ „  . 


or  by  the  Rodrigues  formula: 


Unix)  = 

dx" 

Without  loss  of  generality,  we  will  specialize  to  the  case  t  —  1/4,  and  denote 
fli/4  by  T.  We  now  prove  the  aforementioned 

Observation  1:  T  is  closed  on  P#. 

Proof:  We  will  show  that  THniPn  for 

VtT(  THn)(x)  =  /“  r  ~*P H„(x)dx 

-  f  '  r  '':r2,'(-l)"—  tr-*'  u!\ 
dx" 


I’brc  f> 
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I)""1 


dn~x 

dxn~x 


( e~x2)dx 


=  ViT  2y(Tf/„_i)(y); 


Solving  this  recursion  relation,  using  TH o  =  1,  we  have 

(Tf/„)(x)  =  2"x" .  ■ 


As  a  result  of  Observation  1,  T  is  an  isomorphism  of  Pn .  The  inverse  of  T  on 
/V  >s  clearly  given  by 

,  N  N  at 

T_1  2  (^/x1)  = 

i  =  0  i  =  0  ‘ 

Our  main  result  is  that  T-1  restricted  to  can  be  represented  by  a  convolu¬ 
tion  with  an  explicit  kernel  D«/x): 

Theorem:  For  fiPs  and  j g  =  Tf,  then 

f  =  DN  *  g 


where 


ZV(x)  =  <?  *' 


IJV/2J 

2 


*  =  0 


(-1)* 
ViT  k!2* 


We  will  give  a  proof  below  using  direct  integration  (as  opposed  to  using  Fourier 
transform  distributions).  Note,  however,  that  D/y(x)  is  not  the  unique  function 
representing  T_1  on  PS-  In  general,  the  kernel  can  be  translated  by  any  function 
which  yields  a  zero  convolution  against  Ps ■  This  includes  all  functions  of  the  form 

e~x2Hn(x),  n> N. 

The  stated  kernel  is  unique  among  the  class  of  functions  of  the  form  e~x~P(x), 
where  P(x)  is  a  polynomial  of  degree  N. 


It  is  interesting  to  compare  the  form  of  Dfj(x)  with  standard  enhancement 
filters.  For  example,  for  N  =  3, 


D3(x)  =  ~re  x~ ( 1  x2)—  —~e  x 

VtT  Vtt 


1  d2  1 

2  dx2 


)■ 


Thus 


Oj,  *  g 


\_d]_ 
2  dx2 


Tg. 


which  is  a  not  uncommon  high  emphasis  (liter  (see.  e  g  ,  the  papers  by  H.  Mach 
in  [17]  and  [18].  In  Figure  1,  we  display  plots  of  D\  for  several  values  of  A' 
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Figure  Id 

Figure  1.  The  deblurring  kernels  Ds,  for  N  =  5,9,  and  13  for  Figures  la,  lb,  and 
lc  respectively.  Figures  la,  lb,  and  lc  are  drawn  on  the  same  vertical  scales.  Fig¬ 
ure  Id  shows  the  deblurring  kernel  Z)5  (the  same  as  figure  la)  with  a  different  verti¬ 
cal  scale  ^o  emphasize  the  structure.  The  corresponding  blurring  kernel  is 
(l/ViT)*  ,  so  that  one  standard  deviation  cr  is  equal  to  x=  ±  VT/2. 


The  proof  of  the  theorem  depends  on  several  lemmas. 
Lemma  1: 


Lemma  2: 


An 


f“  —~e~xZxndx 
Vtt 


0, 

n ! 

, 2n(n /2) ! ’ 


n  odd, 
n  even 


C  2k, 2p 


— x‘Hzk(  x)x2pdx 

Vtt 


o. 

p  <k 

_ (Ini' _ 

2 2r~Zkip  ~k  i ’ 

r-K 

Proofs:  Lemma  1  and  Lemma  2  may  be  found  in  [19],  Section  -  16.  Fxampie  1 
Lemma  3: 

For  n  Sr- k. 


<lk  --  f  D\(i)xtdi 


k ' 


: *u 


k  Odd 
k  v  \  c  r. 


Proof:  \  •or  k  odd.  and  using  the  definition  of  \m-  eSivu  r  />«.  >  n‘ 

an  odd  integrable  function,  and  so  integrates  to  zero  For  k  2;-. 
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IA"2J  ,_U1 


f  DN(x)xkdx  =  f  e~x  2  ~  V-  H2i(x)x2pdx 


1-0  i-0  i!2'  22p~2i(j>  —  i)! 

2pp\  £0i'Xp-i)'.  2 Pp\  2 2pP\ 

Proof  of  the  Theorem:  From  the  formula  for  THn  computed  in  Observation  1,  it 
suffices  to  show  that  Ds  *  (2"x")  =  H„(x),n  2:  N.  We  have 


(D*  *  2"x")(>0  =  /"  2 nDS’(x)(y-x)ndx  =  /  2"D*(x)2(-l) 

—  oc  —  oc  - 

1  =  0 


‘My"  k*kc 


n  1 


ZD^.?n . /_nml2wlL  .-2m 


__  v  z  n  ■  (  i  \k.  n-k  j  =  «  t  V  - * - 12 - = _ (  1  \^rnJ  •  v 

*?o  k\(n  —  k)\  y  *  J?0  (2m)  !(n  —  2m) !  22mrn  !  > 


= s  •.•L-Lvr~2m>,"~2m  =  Hn(y)-  ■ 

m  =  om'n  2m). 

The  theorem  could  have  been  proved  using  the  convolution  theorem  and  by 
computing  the  Fourier  transform  of  D/y(x).  Although  we  have  not  taken  this 
approach,  w'e  will  nonetheless  compute  Ds  in  order  to  show  that  the  multiplier  for 
DN  approaches,  pointwise,  the  inverse  of  the  multiplier  for  the  operator  T,  i.e.. 

Observation  2: 

D.v(ai)  -  e“‘/4  pointwise  as  A’- 


Proof: 


»  lV/2J  /_]■)*  ; 

Ds(u)  =  2  1  r-‘  k  X‘Hlk(x)](u), 

k%  VtT  12* 


where  7  stands  for  the  Fourier  transform  operator.  Now, 


7  [e'*:H:t(x)}(  oi)  =  7[(-\)2kJ~ 

dx“ 


\e  *‘)](u>)  =  (iuO^VtT  e 


Thus 


D.v(uj)  =  Y 


l^2\  (_n* 


0  \  TT  A  ’ ' 


•(-  1  )*cd‘*  v it  e 


-  IV  -I 

-U r  ■> 


i  ’  I 

*-<>  *  l  - 


Hence 


1  T  \  .  -  ,  til  *  ** 

hm  -  c  e  -  r 
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It  is  interesting  to  observe  that  the  kernel  D n  is  a  multiplier  where  Fourier 
Transform  is  the  Taylor  series  approximation  to  the  function  e“2/2,  multiplied  by 
the  window  e““2/4. 

Also  as  a  consequence  of  observation  2,  we  see  that  Df/(x)  does  not  converge 
pointwise  to  any  function  as  N  -  since  otherwise  the  Fourier  transform  of  that 
function  would  be  «w2/4,  which  is  impossible.  Dn(x)  does  converge  in  L2(e~x  dx), 
but  that  does  not  imply  pointwise  convergence  to  any  function.  We  accordingly 
have  stable  deblurring  when  using  the  kernels  D^ix),  where  stability  is  measured  in 
terms  of  deviation  from  a  polynomial  of  degree  N ,  and  the  L2(e~x  dx)  norm  is  used 
as  the  metric. 


5.  Higher  Dimensions 

The  Gaussian  blur  operator  is  given  by 

Tj fix)  =  Sm^'nl2e'(x~y)2Ay)dy. 

Due  to  the  separability  of  the  kernel  and  Fubini’s  theorem,  T  can  be  decomposed 
into  n  iterated  blurrings: 

T  =  Tj  o  T2  0  •  •  •  °  T„ 

(Tif)(x)  =  /  ——re~ix‘~y,)2f(xu  ■  ■  ■  ,yit  •  •  ■  ,xn)dy( 

VlT 

Consider  a  polynomial  in  IR": 

/(*)  =  E  aaxa 

|a|s  N 


a  =  (a1(a2,  •  •  •  ,a„),  a,€Z,  a^O,  |a|  =  =  '  '  •  *“"• 

For  fixed  n,  the  function  of  one  real  variable 

gi>'i)  =  fix i,  ■■  ■  ■  •  •  ,xn) 

is  a  polynomial  of  degree  no  greater  than  A\  so 

Ds’  *  (Tg)  =  g 

where  T  is  the  standard  one  dimensional  blurring  operator  (2.1).  Combining  we 
find  that 

fix)  =  SmnDsiyi)Ds{\2)  ■  ■  ■  Dsiyn)iTf)(x-y)dy 

for  any  multivariate  polynomial  fix).  Thus  deblurring  of  blurred  polynomials  of 
degree  N  can  be  accomplished  by  convolution  against  the  kernel 

D\-(x)  =  D\(  I  j  )D\(  i ;  )  D,\  (  ) 


Thus  the  situation  in  higher  dimensions  is  similar  to  the  one  dimensional  case 
The  dcblurring  convolution  kernel  is  separable,  and  will  be  of  the  form  e"‘‘P  tv1, 
where  P(x)  is  a  polynomial  of  degree  nX  in  vCIK"  Figure  2  shows  a  plot  of  D \  tor 
n  —  2,  X  =  3.  Using  a  separable  kernel  tor  dcblurring  has  computation.-.!  advan¬ 
tages,  but  can  lead  to  artifacts  in  diagonal  directions  when  applied  to  certain  images 


P«RC  I  I 


DEBLURRING  GAUSSIAN  BLUR 

We  will  see  these  effects  in  the  next  section. 

6.  Experiments  In  Deblurring 

We  experimented  with  the  deblurring  kernels  Dn  using  both  modest  and  large 
amounts  of  blur.  An  original  image  (Figure  2a)  has  eight  bits  of  grayscale  informa¬ 
tion  digitized  on  a  480  by  S12  grid.  The  image  is  regarded  as  a  piecewise  constant 
function  of  two  continuous  real  variables,  and  blurred  by  convolution  against  the 
Gaussian  {\l-n)e~x  ~y  ,  where  for  Figure  2b  the  interpixel  spacing  is  taken  to  be 
h  =  .15  (corresponding  to  a  standard  deviation  a  =  4.7  pixels).  Images  (2c),  (2d), 
and  (2e)  show  the  results  of  applying  the  deblurring  kernels  D5,  D 9,  and  D 13 
respectively.  The  diagonal  artifacts  arise  due  to  the  use  of  separable  kernels,  and 
become  more  pronounced  for  higher  N.  The  computations  were  done  with  12-bit 
fixed  point  arithmetic  on  a  VICOM  image  processing  computer.  The  results  by 
using  floating  point  arithmetic  on  a  general  purpose  computer  (a  VAX)  were  essen¬ 
tially  identical,  although  the  effects  of  the  diagonal  artifacts  is  very  slightly  reduced. 

For  comparison,  we  display  in  Figure  3  the  results  of  convolving  the  blurred 
image  in  Figure  2b  by  a  kernel  whose  Fou/ier  transform  is  a  truncated  version 
of  e“  .  The  kernel  size  was  arbitrarily  limited  to  100  by  100  pixels,  and  p  was 
chosen  small  enough  so  that  the  kernel  elements  were  considerably  smaller  on  the 
periphery  of  the  kernel.  If  p  were  chosen  too  large,  the  magnitude  of  the  kernel 
elements  decay  extremely  slowly.  Clearly,  deblurring  by  the  kernel  DN  is  far  supe¬ 
rior,  due  to  the  elimination  of  the  sharp  cut-off  in  its  spectral  properties.  The  diffi¬ 
culty  with  a  sharp  cut-off  is  that  it  leads  to  “ringing”  of  the  spatial  kernel,  and  evi¬ 
dence  of  ringing  can  be  seen  in  the  example  of  deblurring. 

In  Figure  4,  we  show  the  blurred  image  using  an  interpixel  distance  of  h  =0.1 
(corresponding  to  a  standard  deviation  a  =  7.1  pixels).  The  blurring  kernel,  to 
beyond  three  standard  deviations,  is  a  roughly  50  by  50  mask.  Figure  4b  shows  this 
image  deblurred  using  the  kernel  Dn 

The  experimental  results  show  that  deblurring  certainly  improves  blurred 
images,  even  with  the  given  inherent  lack  of  stability.  Better  results  may  be  obtain¬ 
able  with  nonlinear  or  stochastic  techniques,  although  our  experiments  certainly 
demonstrate  the  human  visual  system’s  sensitivity  to  visual  quality  of  deblurrings. 


DEBLURRING  GAUSSIAN  BLUR 


Hummel,  Kimia,  and  Zucker 


Figure  2e 


Figure  2.  An  original  image  (2a)  digitized  to  4M>  by  M2  pixels  In  (2b),  the  ima, 
has  been  blurred  by  a  Gaussian  (1  vit'  ‘  ,  with  the  interpixcl  distance  h  0  1 
Figures  (2c),  (2d),  and  (2c)  show  Figure  (2)’)  restored  using  the  deblurring  kerne 
!)<,  D  si,  and  respectiseh 
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I- inure  3.  T  he  blurred  image  in  Figure  (2b)  dcblurrcd  using  a  kernel  whose  Fourier 
Iransform  is  a  truncated  version  of  the  dcblurring  multiplier.  1'hc  result  is  a  pseu- 
doinverse  under  the  blurring  operator.  Difficulties  arise  because  the  deMurring  ker¬ 
nel  "rings." 
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Figure  4b 

Figure  4.  The  original  image  blurred  by  a  Gaussian  with  interpixel  distance 
/i=0.10,  and  the  result  of  dcblurring  using  the  kernel  D  jj. 
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